pkgs <-c("tidyverse", "tidylog", "rTPC", "nls.multstart", "broom", "RColorBrewer", "viridis", "minpack.lm", "patchwork")## minpack.lm needed if using nlsLM()if(length(setdiff(pkgs, rownames(installed.packages()))) >0){install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T) }invisible(lapply(pkgs, library,character.only = T))#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──#> ✔ dplyr 1.1.2 ✔ readr 2.1.4#> ✔ forcats 1.0.0 ✔ stringr 1.5.0#> ✔ ggplot2 3.4.2 ✔ tibble 3.2.1#> ✔ lubridate 1.9.2 ✔ tidyr 1.3.0#> ✔ purrr 1.0.1 #> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──#> ✖ dplyr::filter() masks stats::filter()#> ✖ dplyr::lag() masks stats::lag()#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors#> #> Attaching package: 'tidylog'#> #> #> The following objects are masked from 'package:dplyr':#> #> add_count, add_tally, anti_join, count, distinct, distinct_all,#> distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,#> full_join, group_by, group_by_all, group_by_at, group_by_if,#> inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,#> relocate, rename, rename_all, rename_at, rename_if, rename_with,#> right_join, sample_frac, sample_n, select, select_all, select_at,#> select_if, semi_join, slice, slice_head, slice_max, slice_min,#> slice_sample, slice_tail, summarise, summarise_all, summarise_at,#> summarise_if, summarize, summarize_all, summarize_at, summarize_if,#> tally, top_frac, top_n, transmute, transmute_all, transmute_at,#> transmute_if, ungroup#> #> #> The following objects are masked from 'package:tidyr':#> #> drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,#> spread, uncount#> #> #> The following object is masked from 'package:stats':#> #> filter#> #> #> Loading required package: viridisLite# devtools::install_github("seananderson/ggsidekick") ## not on CRAN library(ggsidekick)theme_set(theme_sleek())# Check the temperature model script! This is the order based on mean temperature, which makes sense for the sharpe schoolfield plot, and therefore we might keep it across plotsorder <-c("SI_HA", "BT", "TH", "SI_EK", "FM", "VN", "JM", "MU", "FB", "BS", "HO", "RA")nareas <-length(unique(order))colors <-colorRampPalette(brewer.pal(name ="RdYlBu", n =10))(nareas)# Load functionshome <- here::here()fxn <-list.files(paste0(home, "/R/functions"))invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))
Read and trim data
Code
d <-#read_csv(paste0(home, "/data/for_analysis/dat.csv")) |> read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for_analysis/dat.csv") |>filter(age_ring =="Y") # use only length-at-age by filtering on age_ring#> Rows: 364546 Columns: 11#> ── Column specification ────────────────────────────────────────────────────────#> Delimiter: ","#> chr (5): age_ring, area, gear, ID, sex#> dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length#> #> ℹ Use `spec()` to retrieve the full column specification for this data.#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#> filter: removed 28,878 rows (8%), 335,668 rows remaining# Sample size by area and cohortns <- d |>group_by(cohort, area) |>summarise(n =n())#> group_by: 2 grouping variables (cohort, area)#> summarise: now 505 rows and 3 columns, one group variable remaining (cohort)# Minimum number of observations per area and cohortd$area_cohort <-as.factor(paste(d$area, d$cohort))d <- d |>group_by(area_cohort) |>filter(n() >100)#> group_by: one grouping variable (area_cohort)#> filter (grouped): removed 2,844 rows (1%), 332,824 rows remaining# Minimum number of observations per area, cohort, and aged$area_cohort_age <-as.factor(paste(d$area, d$cohort, d$age_bc))d <- d |>group_by(area_cohort_age) |>filter(n() >10)#> group_by: one grouping variable (area_cohort_age)#> filter (grouped): removed 3,974 rows (1%), 328,850 rows remaining# Minimum number of cohorts in a given areacnt <- d |>group_by(area) |>summarise(n =n_distinct(cohort)) |>filter(n >=10)#> group_by: one grouping variable (area)#> summarise: now 12 rows and 2 columns, ungrouped#> filter: no rows removedd <- d[d$area %in% cnt$area, ]# Plot cleaned dataggplot(d, aes(age_bc, length_mm)) +geom_jitter(size =0.1, height =0, alpha =0.1) +scale_x_continuous(breaks =seq(20)) +theme(axis.text.x =element_text(angle =0)) +theme(axis.text =element_text(size =12), axis.title =element_text(size =15)) +labs(x ="Age", y ="Length (mm)") +guides(color ="none") +facet_wrap(~area, scale ="free_x")
Code
# Longitude and latitude attributes for each areaarea <-c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH", "VN")nareas <-length(area)lat <-c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1, 57.5)lon <-c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9, 16.9)area_attr <-data.frame(cbind(area = area, lat = lat, lon = lon)) |>mutate_at(c("lat","lon"), as.numeric)#> mutate_at: converted 'lat' from character to double (0 new NA)#> converted 'lon' from character to double (0 new NA)
Fit VBGE models
Code
# Get individual growth parameters (functions: VBGF/Gompertz,nls_out,fit_nls)IVBG <- d |>group_by(ID) |>summarize(k =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$k,k_se =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$k_se,linf =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$linf,linf_se =nls_out(fit_nls(length_mm, age_bc, min_nage =4, model ="VBGF"))$linf_se)#> group_by: one grouping variable (ID)#> summarize: now 91,422 rows and 5 columns, ungrouped
Inspect predictions
Code
IVBG <- IVBG |>drop_na(k) # The NAs are because the model didn't converge or because they were below the threshold age#> drop_na: removed 50,514 rows (55%), 40,908 rows remainingIVBG <- IVBG |>mutate(k_lwr = k -1.96*k_se,k_upr = k +1.96*k_se,linf_lwr = linf -1.96*linf_se,linf_upr = linf +1.96*linf_se,row_id =row_number())#> mutate: new variable 'k_lwr' (double) with 40,864 unique values and 0% NA#> new variable 'k_upr' (double) with 40,864 unique values and 0% NA#> new variable 'linf_lwr' (double) with 40,864 unique values and 0% NA#> new variable 'linf_upr' (double) with 40,864 unique values and 0% NA#> new variable 'row_id' (integer) with 40,908 unique values and 0% NA# Plot all K'sIVBG |>#filter(row_id < 5000) |>ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr)) +geom_point(alpha =0.5) +geom_errorbar(alpha =0.5) +NULL
# We can also consider removing the upper 95% quantile of standard errors (not quantile of K directly)IVBG |> tidylog::mutate(keep =ifelse(k >quantile(k_se, probs =0.95), "N", "Y")) |>#filter(row_id < 10000) |>ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +geom_point(alpha =0.5) +facet_wrap(~keep, ncol =1) +geom_errorbar(alpha =0.5) +NULL#> mutate: new variable 'keep' (character) with 2 unique values and 0% NA
Code
# Add back cohort and area variablesIVBG <- IVBG |>left_join(d |>select(ID, area_cohort)) |>separate(area_cohort, into =c("area", "cohort"), remove =TRUE, sep =" ") |>mutate(cohort =as.numeric(cohort))#> Adding missing grouping variables: `area_cohort_age`#> select: dropped 10 variables (length_mm, age_bc, age_catch, age_ring, area, …)#> Joining with `by = join_by(ID)`#> left_join: added 2 columns (area_cohort_age, area_cohort)#> > rows only in x 0#> > rows only in y (115,512)#> > matched rows 213,338 (includes duplicates)#> > =========#> > rows total 213,338#> mutate: converted 'cohort' from character to double (0 new NA)# Summarise and save for sample sizesample_size <- IVBG |>group_by(area) |>summarise(n_cohort =length(unique(cohort)),min_cohort =min(cohort),max_cohort =min(cohort),n_individuals =length(unique(ID)),n_data_points =n())#> group_by: one grouping variable (area)#> summarise: now 12 rows and 6 columns, ungroupedwrite_csv(sample_size, paste0(home, "/output/sample_size.csv"))# Compare how the means and quantiles differ depending on this filteringIVBG_filter <- IVBG |> tidylog::drop_na(k_se) |> tidylog::filter(k_se <quantile(k_se, probs =0.95))#> drop_na: no rows removed#> filter: removed 10,670 rows (5%), 202,668 rows remaining# Summarize growth coefficients by cohort and areaVBG <- IVBG |>group_by(cohort, area) |>summarize(k_mean =mean(k, na.rm = T),k_median =quantile(k, prob =0.5, na.rm = T),linf_median =quantile(linf, prob =0.5, na.rm = T),k_lwr =quantile(k, prob =0.05, na.rm = T),k_upr =quantile(k, prob =0.95, na.rm = T)) |>ungroup()#> group_by: 2 grouping variables (cohort, area)#> summarize: now 382 rows and 7 columns, one group variable remaining (cohort)#> ungroup: no grouping variablesVBG_filter <- IVBG_filter |>group_by(cohort, area) |>summarize(k_mean =mean(k, na.rm = T),k_median =quantile(k, prob =0.5, na.rm = T),k_lwr =quantile(k, prob =0.05, na.rm = T),k_upr =quantile(k, prob =0.95, na.rm = T)) |>ungroup()#> group_by: 2 grouping variables (cohort, area)#> summarize: now 382 rows and 6 columns, one group variable remaining (cohort)#> ungroup: no grouping variablesggplot() +geom_ribbon(data = VBG, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr, fill ="All k's"), alpha =0.4) +geom_ribbon(data = VBG_filter, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr, fill ="Filtered"), alpha =0.4) +geom_line(data = VBG, aes(cohort, k_median, color ="All k's")) +geom_line(data = VBG_filter, aes(cohort, k_median, color ="Filtered")) +guides(fill =FALSE) +facet_wrap(~area)#> Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as#> of ggplot2 3.3.4.
Code
ggplot() +geom_line(data = VBG, aes(cohort, k_median, color ="All k's")) +geom_line(data = VBG_filter, aes(cohort, k_median, color ="Filtered")) +guides(fill =FALSE) +facet_wrap(~area)
Code
# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.
Add GAM-predicted temperature to growth models
Code
pred_temp <-read_csv(paste0(home, "/output/gam_predicted_temps.csv")) |>rename(cohort = year)#> Rows: 566 Columns: 4#> ── Column specification ────────────────────────────────────────────────────────#> Delimiter: ","#> chr (2): area, type#> dbl (2): year, temp#> #> ℹ Use `spec()` to retrieve the full column specification for this data.#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#> rename: renamed one variable (cohort)VBG <- VBG |>left_join(pred_temp, by =c("area", "cohort"))#> left_join: added 2 columns (temp, type)#> > rows only in x 0#> > rows only in y (184)#> > matched rows 382#> > =====#> > rows total 382# save data for map-plotcohort_sample_size <- IVBG |>group_by(area, cohort) |>summarise(n =n()) # individuals, not samples!#> group_by: 2 grouping variables (area, cohort)#> summarise: now 382 rows and 3 columns, one group variable remaining (area)VBG <-left_join(VBG, cohort_sample_size, by =c("cohort", "area"))#> left_join: added one column (n)#> > rows only in x 0#> > rows only in y ( 0)#> > matched rows 382#> > =====#> > rows total 382write_csv(VBG, paste0(home, "/output/vbg.csv"))
Plot VBGE fits
Code
ids <- d |>group_by(area) |>sample_n(50)#> group_by: one grouping variable (area)#> sample_n (grouped): removed 328,250 rows (>99%), 600 rows remainingdat <- IVBG %>%filter(ID %in% ids$ID)#> filter: removed 211,573 rows (99%), 1,765 rows remainingd |>ungroup() |>filter(ID %in% ids$ID) |>left_join(dat) |> tidylog::mutate(length_mm_pred = linf*(1-exp(-k*age_bc))) |>ggplot(aes(age_bc, length_mm, group = ID, color = ID)) +geom_jitter(width =0.2, height =0, alpha =0.5, size =0.6) +geom_line(aes(age_bc, length_mm_pred, group = ID, color = ID), inherit.aes =FALSE, alpha =0.5) +guides(color ="none") +scale_color_viridis(discrete =TRUE, name ="Area", option ="cividis") +labs(x ="Age (yrs)", y ="Length (mm)") +facet_wrap(~factor(area, order))#> ungroup: no grouping variables#> filter: removed 326,388 rows (99%), 2,462 rows remaining#> Joining with `by = join_by(area, cohort, ID, area_cohort_age)`#> left_join: added 9 columns (k, k_se, linf, linf_se, k_lwr, …)#> > rows only in x 697#> > rows only in y ( 0)#> > matched rows 1,765#> > =======#> > rows total 2,462#> mutate: new variable 'length_mm_pred' (double) with 1,766 unique values and 28%#> NA#> Warning: Removed 697 rows containing missing values (`geom_line()`).
Code
ggsave(paste0(home, "/figures/vb_fits.pdf" ), width =17, height =17, units ="cm")#> Warning: Removed 697 rows containing missing values (`geom_line()`).k <- IVBG |>ggplot(aes(k, color =factor(area, order))) +geom_density(alpha =0.4, fill =NA) +scale_color_viridis(discrete =TRUE, name ="Area", direction =-1) +coord_cartesian(expand =0) +guides(color ="none") +theme(legend.position =c(0.8, 0.8))linf <- IVBG |>ggplot(aes(linf, color =factor(area, order))) +geom_density(alpha =0.4, fill =NA) +scale_color_viridis(discrete =TRUE, name ="Area", direction =-1) +coord_cartesian(expand =0, xlim =c(0, 2000)) +theme(legend.position =c(0.9, 0.5))k / linf
Code
ggsave(paste0(home, "/figures/vb_pars.pdf" ), width =17, height =17, units ="cm")
Fit Sharpe-Schoolfield model to K
By area
Code
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat <- VBG |>select(k_median, temp) |>rename(rate = k_median)#> select: dropped 8 variables (cohort, area, k_mean, linf_median, k_lwr, …)#> rename: renamed one variable (rate)lower <-get_lower_lims(dat$temp, dat$rate, model_name = model)upper <-get_upper_lims(dat$temp, dat$rate, model_name = model)start <-get_start_vals(dat$temp, dat$rate, model_name = model)# Sharpe-Schoolfield model fit to data for each areapreds <-NULLfor(a inunique(VBG$area)) {# Get data dat <- VBG |>filter(area == a) |>select(k_median, temp, area) |>rename(rate = k_median)# Fit model fit <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new data new_data <-data.frame(temp =seq(min(dat$temp), max(dat$temp), length.out =100)) pred <-augment(fit, newdata = new_data) |>mutate(area = a)# Add to general data frame preds <-data.frame(rbind(preds, pred))}#> filter: removed 319 rows (84%), 63 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 328 rows (86%), 54 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 341 rows (89%), 41 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 347 rows (91%), 35 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 334 rows (87%), 48 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 344 rows (90%), 38 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 363 rows (95%), 19 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 348 rows (91%), 34 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 364 rows (95%), 18 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 367 rows (96%), 15 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 370 rows (97%), 12 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 377 rows (99%), 5 rows remaining#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA
All areas pooled
Code
# One sharpe schoolfield fitted to all datafit_all <-nls_multstart( k_median ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = VBG,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new datanew_data_all <-data.frame(temp =seq(min(VBG$temp),max(VBG$temp),length.out =100))pred_all <-augment(fit_all, newdata = new_data_all) |>mutate(area ="all")#> mutate: new variable 'area' (character) with one unique value and 0% NA# Add t_optkb <-8.62e-05data.frame(par =names(coef(fit_all)),est =coef(fit_all)) |>pivot_wider(names_from = par, values_from = est) |>summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))#> pivot_wider: reorganized (par, est) into (r_tref, e, eh, th) [was 4x2, now 1x4]#> summarise: now one row and one column, ungrouped#> # A tibble: 1 × 1#> t_opt#> <dbl>#> 1 19.6
Plot Sharpe Schoolfield fits
Code
# Plot growth coefficients by year and area against mean SSTpreds |>ggplot(aes(temp, .fitted, color =factor(area, order))) +geom_point(data = VBG, aes(temp, k_median, color =factor(area, order)), size =0.6) +geom_line(aes(temp, .fitted, group =factor(area)), linewidth =1) +geom_line(data = pred_all, aes(temp, .fitted),linewidth =1, inherit.aes =FALSE, linetype =2) +#scale_color_manual(values = colors, name = "Area") +scale_color_viridis(name ="Area", direction =-1, discrete =TRUE) +theme(plot.title =element_text(size =15, face ="bold")) +theme(axis.text.x =element_text(angle =0)) +theme(axis.text =element_text(size =12), axis.title =element_text(size =15)) +guides(color =guide_legend(override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 1)) +labs(x ="Temperature (C)", y ="Median von Bertalanffy growth coefficient, k") +NULL
Code
ggsave(paste0(home, "/figures/sharpe_school.pdf" ), width =17, height =17, units ="cm")
Test plotting Linf aswell…
By area
Code
model <-'sharpeschoolhigh_1981'# Get starting values on full dataset for Sharpe-Schoolfield modeldat <- VBG |>select(linf_median, temp) |>rename(rate = linf_median)#> select: dropped 8 variables (cohort, area, k_mean, k_median, k_lwr, …)#> rename: renamed one variable (rate)lower <-get_lower_lims(dat$temp, dat$rate, model_name = model)upper <-get_upper_lims(dat$temp, dat$rate, model_name = model)start <-get_start_vals(dat$temp, dat$rate, model_name = model)# Sharpe-Schoolfield model fit to data for each areapreds <-NULLfor(a inunique(VBG$area)) {# Get data dat <- VBG |>filter(area == a) |>select(linf_median, temp, area) |>rename(rate = linf_median)# Fit model fit <-nls_multstart( rate ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = dat,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new data new_data <-data.frame(temp =seq(min(dat$temp), max(dat$temp), length.out =100)) pred <-augment(fit, newdata = new_data) |>mutate(area = a)# Add to general data frame preds <-data.frame(rbind(preds, pred))}#> filter: removed 319 rows (84%), 63 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 328 rows (86%), 54 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 341 rows (89%), 41 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 347 rows (91%), 35 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 334 rows (87%), 48 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 344 rows (90%), 38 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 363 rows (95%), 19 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 348 rows (91%), 34 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 364 rows (95%), 18 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 367 rows (96%), 15 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 370 rows (97%), 12 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA#> filter: removed 377 rows (99%), 5 rows remaining#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)#> rename: renamed one variable (rate)#> mutate: new variable 'area' (character) with one unique value and 0% NA
All areas pooled
Code
# One sharpe schoolfield fitted to all datafit_all <-nls_multstart( linf_median ~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref =8),data = VBG,iter =c(3, 3, 3, 3),start_lower = start*0.5,start_upper = start*2,lower = lower,upper = upper,supp_errors ='Y' )# Make predictions on new datanew_data_all <-data.frame(temp =seq(min(VBG$temp),max(VBG$temp),length.out =100))pred_all <-augment(fit_all, newdata = new_data_all) |>mutate(area ="all")#> mutate: new variable 'area' (character) with one unique value and 0% NA# Add t_opt# kb <- 8.62e-05# data.frame(par = names(coef(fit_all)),# est = coef(fit_all)) |> # pivot_wider(names_from = par, values_from = est) |> # summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))
Plot Sharpe Schoolfield fits
Code
# Nice color palette, order based on average temperature!# Plot growth coefficients by year and area against mean SSTpreds |>ggplot(aes(temp, .fitted, color =factor(area, order))) +geom_point(data = VBG, aes(temp, linf_median, color =factor(area, order)), size =0.6) +geom_line(aes(temp, .fitted, group =factor(area)), linewidth =1) +geom_line(data = pred_all, aes(temp, .fitted),linewidth =1, inherit.aes =FALSE, linetype =2) +# geom_smooth(aes(temp, .fitted, color = factor(area, order)), # method = "gam", formula = y ~ s(x, k = 4)) +# geom_smooth(aes(temp, .fitted), inherit.aes = FALSE,# method = "gam", formula = y ~ s(x, k = 4), linetype = 2) +#scale_color_manual(values = colors, name = "Area") +scale_color_viridis(name ="Area", direction =-1, discrete =TRUE) +theme(plot.title =element_text(size =15, face ="bold")) +theme(axis.text.x =element_text(angle =0)) +theme(axis.text =element_text(size =12), axis.title =element_text(size =15)) +guides(color =guide_legend(override.aes =list(size =1))) +scale_x_continuous(breaks =seq(-5, 20, 1)) +labs(x ="Temperature (C)", y ="Asymptotic size (mm)") +NULL
